Adaptive logarithmic discretization for numerical renormalization group 

methods 



Rok Zitko a 

" Institute for Theoretical Physics, University of Gottingen, 

Friedrich- Hund- Platz 1, D-37077 Gottingen, Germany 
J. Stefan Institute, Jamova 39, SI- 1000 Ljubljana, Slovenia 



Abstract 

The problem of the logarithmic discretization of an arbitrary positive function (such as the density of states) 
is studied in general terms. Logarithmic discretization has arbitrary high resolution around some chosen 
point (such as Fermi level) and it finds application, for example, in the numerical renormalization group 
(NRG) approach to quantum impurity problems (Kondo model), where the continuum of the conduction 
band states needs to be reduced to a finite number of levels with good sampling near the Fermi level. 
The discretization schemes under discussion are required to reproduce the original function after averaging 
over different interleaved discretization meshes, thus systematic deviations which appear in the conventional 
logarithmic discretization are eliminated. An improved scheme is proposed in which the discretization-mesh 
points themselves are determined in an adaptive way; they are denser in the regions where the function 
has higher values. Such schemes help in reducing the residual numeric artefacts in NRG calculations in 
situations where the density of states approaches zero over extended intervals. A reference implementation 
of the solver for the differential equations which determine the full set of discretization coefficients is also 
described. 
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1. Introduction 

Discretization of the continuum of conduction-band electron states to a finite set of levels is a com- 
mon approximation procedure in many practical problems in computational condensed-matter physics. In 
band structure calculations, for example, the discretization is performed in the reciprocal space where a 
finite number of judiciously chosen (symmetry-adapted) crystal- momentum points are chosen to sample 
contributions from various regions of the full Brillouin zone [![. Various real-space approaches based on 
finite-difference or finite-element formulations are also possible 0. Finally, in the numerical renormaliza- 
tion group (NRG) methods for solving quantum impurity problems the quantity that is discretized is the 
density of states, thus the discretization is effectively performed in the energy space [3,0,0,0]. In order to 
obtain good sampling of states near the Fermi level, the discretization mesh consists of a geometric series of 
points, so that very high resolution is achieved in the vicinity of the Fermi level [3|]; this choice is motivated 
by the nature of the quantum impurity problems (described by models such as the Kondo model or the 
Anderson impurity model), where excitations from different energy scales have a comparable effect on the 
physical properties. NRG allows to calculate dynamic pro perties (spectral functions, dynamic susceptibili- 
ties) of impurity models @, 0, 0, El, El, EH, El El El, El, El, E3, El, El El and, in particular, it may be 



used as an impurity solver in the dynamic mean-field theory (DMFT) approach to the strongly-correlated 



electron systems |21|, |22|, |23|, [24|, |25|, |26| . Logarithmic discretization of the conductance band can also be 
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used in the density-matrix renormalization-group (DMRG) calculations [III, [H, El , in the embedded-cluster 
approximation [30( , and in the artificial Friedel resonance approach 31 1 . 

Recent work has revealed that the conventional discretization schemes used in NRG lead to discretization 
artefacts, which restrain the accuracy of calculations and limit the ultimate resolution that may be achieved 
[20I ] . A new approach to the discretization based on requiring exact reproduction of the conduction-band 
density of states (after discretization-mesh averaging) was shown to lead to a notable reduction of the 
artefacts 20]. In this paper the approach is further improved by making the discretization grid itself 
adaptable. This work is organized as follows: in Sec. [2] the discretization problem is formulated in very 
general terms starting from basic principles in order to determine the most general form of the discretization 
equations. In Sec. 02 the origin of the residual artefacts in the fixed-discretization-mesh approach is studied, 
while Sec. [5] presents a simple way of adapting the discretization grid to the density of states. Finally, in 
Sec. [5] the implementation of the discretization equation solver is described. 

2. Problem statement 

We consider a non-interacting Hamiltonian for the conductance-band electrons describing a set of N 
states indexed by some quantum number k: 

N 



H = J2tk\k)(k\. (1) 



fc=i 

We assume that for all k, Ck G [— D : D], and in the following we use D as the unit of energy. The density 
of states at energy u> is defined as 

p(") = ^Y, 6 ^ ~ e ^ ( 2 ) 

k 

and it is normalized to 1. In the continuum N — > 00 limit, p(u>) is an arbitrary positive function with finite 
support [—1 : 1]. Our goal is to find a discrete representation of the continuum of states that is suitable 
for numeric calculations in quantum impurity physics and which reproduces the density of states p(u>) as 
accurately as possible. The approach is clearly applicable to an arbitrary positive function (in the single- 
impurity Anderson model, for example, the function of interest is the hybridisation function T), but to be 
specific we may think of it as a density of states. 

In the following we focus on positive energies u>; for negative uj the procedure is fully equivalent. We 
discretize the interval [0 : 1] by defining the discretization- mesh points ej , j € N, such that 

ei = l, 6j > Cj+i, lim Ej — 0. (3) 

j—tOO 

In order to achieve high resolution near lo = (i.e. in the vicinity of the Fermi level), we furthermore require 
that 6j behave asymptotically as 

ej ~A-'\ (4) 

where A > 1 is a constant number known as the discretization parameter. We denote each discretization 
interval [ej+i ■ ej] as Ij. Within each interval Ij we choose a representative energy £j. Approximating 
£j ej, we see that the relative "resolution" 

Ag „ Ej-E j+1 _ A-l 

E ~ Ej + E j+1 ~ A + 1 [ ' 

is approximately the same on all energy scales; from this property stems the name "logarithmic discretiza- 
tion". We furthermore notice that the resolution is improved as A is reduced toward A = 1, which corre- 
sponds to returning to the continuum limit. At given A greater than 1, however, the continuum spectral 
density is represented by a set of delta peaks, 

p{u) — WjS(uj — £j), for lj > 0, (6) 

j 



where weights Wj need to be chosen so that the normalization is preserved: 

p{oj)duj = W. (7) 



Here we have introduced the total spectral weight for positive frequencies, W. The prescription for determin- 
ing the discretization- mesh points Cj , representative energies £j and weights Wj (and their negative-frequency 
analogues e J , £j , and wj ) is what we refer to as a "discretization scheme" . Several such procedures are 
known from the literature on the numerical renormalization group @, S m, M M M M, Ha M, EqI . 



A commonly used technique to improve the accuracy of numerical calculations is to perform computations 
for several different discretizations of the continuum and average the results. In the context of the numerical 
renormalization group, this is known as the "z-averaging" or the "interleaved method" [^, 3^, IB]. In this 



approach one introduces a continuous parameter z 6 [0 : 1] which characterizes different discretizations 
{cj, £j, w*j}. Parameter z is sometimes called the "twist parameter" since it can be related to twisting the 
boundary conditions of the wave-functions of the conduction-band electrons [l6| . The z-averaging procedure 
is meaningful when the boundary conditions 

Vz el = 1, S° = 1, (8) 

and continuity constraints 

Vj e)=e° j+l , Vj £]=£° j+1 , (9) 

are satisfied. We also require that e| and £* be monotonously decreasing functions of z, but they need not 
be continuous. Furthermore, for each value of z, the weights need to be normalized: 

oo 

5>< = W. (10) 

3=1 

The spectral density is then determined as the average (integral) over z [l^ : 

= [ p(z,w)dz, (11) 



o 



where p(z, us) is given by the generalization of Eq. 



3 

We obtain M 



A(u) = Y; l w*5(u;-q)dz = -^L_, (13) 
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where the parameters j and z in the last part of the expression are determined implicitly through the relation 
lu = £j] this choice of j and z is unique due to the requirement that f J is strictly decreasing as a function 
of z. In practical calculations, only a small number N z of values of z is used (often as small as 4 or 2, 
or even a single one) and the integral is approximated using the rectangle method; the resulting function 
then needs to be broadened appropriately to obtain a continuous representation. See, however, Ref. for 
calculations with a large N z and very narrow broadening kernel, which show that overbroadening errors of 
NRG can largely be eliminated. 

To simplify the notation and to provide more insight into the mathematical structure of the problem, 
we introduce continuous indexing of the discretization-grid points as [20| 



x = j + z. 
3 



(14) 



The "grid parameter" x runs from 1 to +00 and the coefficients e|, £J and io| become continuous functions 
e(a;), £{x) and u?(a;), while the continuity constraints, Eq. ([5]), are automatically satisfied. Both e(x) and 
£ (x) are required to be monotonously decreasing and, furthermore, the following boundary conditions need 
to be satisfied: 

e(x) = 1 for x < 2, lim e{x) = 0, (15) 

x — >oo 

£(1) = 1, lim5(i)=0. (16) 

x — >oo 

These equations embody the requirement that as the grid parameter x sweeps the interval [1 : +00), both 
the grid energy e(x) and the representative energy £ (x) describe the totality of the unit interval [0 : 1]. For 
convenience, we define the inverse function R to £ as x = R(ui), i.e. Ro £ = 1. 
Using newly introduced notation, Eq. (|13p becomes 

A ^ = ^W^x^ ithx = R ^ < 17 > 

The normalization condition can be expressed as 

A{cu)duj = / A(iu)=-^dx= / w{x)Ax^W 1 (18) 

Joo dx J 1 

but, in addition, the more general normalization equation 

00 

^w(j+z) = W, (19) 
j=i 

which follows from Eq. (flU)) . must be satisfied for each z E [0 : 1]; in fact, Eq. (fT5| follows trivially from 
Eq. (|19p. We notice that the function e(x) appears neither in Eq. (fTTf nor in Eq. (fTO| : grid coefficients e(x) 
are thus only auxiliary quantities which define the discretization grid without actually explicitly appearing 
in the final result for the density of states. 

With all these preparations we can now finally state the problem as follows: we seek to determine 
functions £{x) and w{x) such that 

w\R(u>)] , , 

-6£[rI,)\/4* =P{u) ' (20) 

and satisfying normalisation, Eq. (|19|) . 

Since, e(x) = 1 for x € [1 : 2] and e(x) ->0asi-> 00, it is easy to see that w(x) defined as 

w(x) = / p(uj)duj (21) 

Je(x+1) 

solves the normalization Eq. (fT5|) . In fact, this is the only general solution of Eq. (|19p : we note, however, that 
e(x) can in principle still be an arbitrary monotonously decreasing function on the interval [2 : 00). Thus w(x) 
must be defined as the integral of the density of states in the discretization integral I(x) = [e(x + 1) : e(x)), 
while we still have the full liberty of choosing the discretization mesh in any convenient way. 

To make contact with the original motivation for introducing the logarithmic discretization of a contin- 
uum, we now focus on discretization meshes e(x) with asymptotic behavior e(x) ~ A 2 ~ x , where we have 
shifted the exponent by 2 for convenience. This asymptotic form is equivalent to requiring 

^1 = -A 2 - \nAC(x,e), (22) 
da; 

where C(x, e) is an arbitrary strictly positive function with a non-zero limit 

A=lim C(x,e), (23) 

4 



The requirement A 7^ is necessary to obtain the desired asymptotic behavior. 

The problem is thus reduced to finding an appropriate function C(x, e). Once C(x, e) is chosen, the full 
solution of the problem is obtained by solving the initial value problem 



8(1) = 1, 

dx p[£(x)] 

In NRG, the discrete levels resulting from the discretization for given twist parameter z are used to write 
the discretized form of the conduction-band Hamiltonian 

00 

H = J2mi)(i\- (25) 

i=l 

One then forms the combination of levels 

I/o)=EM. (26) 

i 

and transforms the Hamiltonian H into a new basis (the first state of which is |/o)), so that the Hamiltonian 
takes the form of a tight-binding chain 0, d, l38l. Iti|: 

oo oo 

H = £»|/n></»| + E *» d/«)(/n+l| + |/n+l)(/n|) • (27) 
n=Q n=0 

This form is known as the "hopping Hamiltonian" or the "Wilson chain" . The coefficients t n decrease 
asymptotically as t n cx A~"/ 2 and only a finite number of chain sites is retained in practical NRG calculations. 
The impurity hybridizes with the conduction band through level |/o) only; this implies that the spectral 
function Af (lu) of the level /o represents the density of states of the conduction band as seen by the impurity. 
The goal is thus to make Af(u) a good representation of the density of states p(to). 

We will discuss possible choices of C(x,e) in Sec. 3] and describe a numerical approach for solving the 
initial value problem in Sec. but we first pause to describe the previously known discretization schemes 
based on fixed discretization mesh to point out the origin of their deficiencies. 



3. Fixed discretization mesh 

The conventional logarithmic discretization used in NRG calculations is based on the fixed discretization 
mesh obtained by setting C(x, e) = 1 in Eq. (|22|) . with solution 

The different discretization schemes then differ in the recipe for calculating £{x). The first schemes for an 
arbitrary density of states [H, [H, HI] used a physically motivated, but otherwise ad-hoc expression 

£( X ) = M *£ )PK ' , (29) 

?H dco 

i.e. an average of energy weighted by p(u>). An improved approach was later introduced [l6| . where 

f ( X ) = ' — . (30) 
/r+i)P( w )/ u(iu 



This scheme has better convergence properties as A — > 1 and it doesn't systematically underestimate the 
density of states at low energies. Neither of these two approaches, however, in general satisfies Eq. (j2"0"l) . 
which leads to artefacts which become apparent in high-resolution NRG calculations, in particular at high 
energies near the band edges ■ More recently, a scheme based on solving Eq. (|24|) for a fixed discretization 
mesh was shown to be very successful in removing the most severe of these artefacts [2(j. In most cases, this 
approach works well, it is robust and the conductance-band density of states can be reproduced accurately 
in numerical calculations. We find, however, that in practical NRG calculations some residual artefacts 
still appear, in particular in situations where the density of states has large and rapid variations. These 
artefacts cannot be fully removed by increasing the number of z- values used in the averaging. The artefacts 
appear especially near energies £j, j € Z. Their origin can be traced back to the systematic errors in NRG 
calculations, which shift the spectral peaks from exact energies to £j+6£j. For given j, the error 6£j is a 

smooth function of parameter z, thus 8£^ Nz — 8£^j n+1 ^ Nz is a small quantity for all n — 1, . . ., N z — 1. Even 
a narrow broadening kernel will then give a smooth final result in these energy ranges and the resolution can 
in principle be systematically improved by decreasing the broadening width while simultaneously increasing 
the total number of z values, N z . To the contrary, the difference in the values of 8£^* and 8£j^ z is 
large, since the two calculations are based on two very different discretization meshes. This is still true 
even for very large N z . Near £j points, thus, there exist a minimum broadening width below which these 
systematic errors become unmasked and manifest themselves as sharp artefacts. This effectively limits the 
highest spectral resolution that one may expect to achieve in NRG calculations. 

It was found that problems of this type become especially pronounced for densities of states which are 
very low (or zero) over considerable energy intervals. Such situation commonly occurs when NRG is used as 
the impurity solver in the dynamical mean-field theory [24|j [2|| [2(| . It should be noted that a representative 
state will be chosen in each interval even if the interval contains very small (or even zero) spectral weight. 
While such states appear with little weights w i: their presence is nevertheless found to be detrimental to 
the overall accuracy of the calculation. The issue is, actually, somewhat subtle: the weights wi only affect 
the combination of states that forms |/o). The terms £i\i)(i\ in regions of small density of states still appear 
in H of Eq. (|27[) and they affect the coefficients £„, t n of the Wilson chain in a way which is detrimental 
for the accuracy of the calculation. This problem is especially severe for very small z, where some of the 
coefficients £ n may become very large (when expressed in units of A~™/ 2 ). 

We now study the severity of this problem on the example of a density of states with a smooth transition 
from lower density at high energies to higher density at lower energies. As a simple model we choose a 
density of states described by 

p(u>) oc I + wtanh [— cl(\lj\ — loq)] , (31) 

where w (0 < < 1) determines the height of the variation, loq the energy where the change occurs and a 
its rapidity. 

In Figs. [1] and Owe compare Af(u>), the spectral function on the first site of the Wilson chain (level /o), 
computed using the NRG, with the model density of states p(co). Calculations were performed with N z = 64 
z values, discretization parameter A = 2, truncation cutoff set at IOlon, broadening parameter r\ = 0.01 and 
patching parameter p = 2.1 (for details on the method see Ref. [1(3]). We focus now on left panels which 
contain results obtained using fixed grid discretization; right panels with adaptive-grid calculation results 
will be discussed in Sec.|H In Fig. [1] the transition is made progressively sharper. At first only the artefacts 
at energies £j are amplified, however for large enough a additional sharp features appear. In Fig. [21 we 
shift the transition point to lower energies, thereby again increasing the range of very low density of states. 
While for ujq — 0.9 the artefacts are very mild, they increase significantly as soon as p(oj) becomes very low 
for large lo. It may be noted that mild artefacts may easily be removed by making the broadening kernel 
wider, but one then needs to renounce on obtaining high-resolution spectral functions. When artefacts are 
severe, they might affect even calculations with very wide broadening kernel. 
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Figure 1: Spectral function of the first site on the Wilson chain for tanh density of states of the conduction band: comparison 
between results using fixed and adaptive discretization grids for a range of parameter a, which controls the steepness of the 
transition and consequently the value of p(u>) in the low-density range. 
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Figure 2: Spectral function of the first site on the Wilson chain for tanh density of states of the conduction band: comparison 
between results using fixed and adaptive discretization grids for a range of parameter luq, which controls the energy of the 
transition from high-density to low-density range. 
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4. Adaptive discretization mesh 



From the considerations detailed in the previous section it follows that a more appropriate choice of the 
discretization grid should be such that the grid were less dense in the regions of low density of states. A 
convenient way to implement this requirement is by demanding that the function C'(x, e) in Eq. (|22p be of 
the form 

C(x,e) = ^-. (32) 

It can be seen that in the limit p(e) — > 0, the corresponding energy regions will not even appear in the 
discretization grid, which is a very desirable property. 

We consider a specialization to f(x) = A, where A is a constant. It is easy to see that A must, in fact, 
be equal to the total weight W = f Q p(uj) duj. We also observe that for a featureless flat band [p(lo) = 1/2], 
this approach reduces to the conventional fixed discretization grid. 

We furthermore note that the "characteristic energy scale at the iVth NRG iteration", ujn, depends on 
the discretization grid (i.e. on function e(x)), thus it becomes dependent on p(lu). This has important 
consequences on the choice of parameters in NRG calculations; in particular, the truncation cutoff and the 
patching parameter p need to be appropriately redefined. A simple choice is to rescale a; at by the ratio 

R= lim i^L (33) 

where €q(x) is the discretization grid function for fixed grid. This choice appears suitable for the tanh test 
functions, but it is not expected to be generally applicable. For strongly varying density of states, the 
patching procedure itself might become problematic and it might be better to use the complete-Fock-space 



approach [15j, [171, [18| , although that technique is not without problems either [20( . The issue of extracting 
spectral information from partial results at different iterations clearly merits further attention. 

Comparing the right panels in Figs. [T] and [5] with the left panels, we see that the adaptive-grid approach 
leads to a considerable improvement; the artefacts essentially disappear. We notice that some oscillations 
appear in the transition region between low and high spectral density, see the inset in Fig. [TJ It is important 
to notice that if the discretization mesh is adaptive, the width of the broadening kernel in the calculation 
of the continuous spectral function should ideally also take the mesh density into account; energy regions 
where mesh is less dense should be broadened more, since there will be less representative energy points. 
At the same time, one should be careful to avoid excessive broadening of spectral functions on the impurity 
levels, since low density in the conductance band translates into sharper features in the impurity spectral 
function. We note that in the extreme case of a gap in the conductance band, there might even appear 
delta peaks in the impurity spectral function. One should thus base the choice of the broadening width on 
physical considerations and on the expected spectral features. 

5. Initial value problem 

Both e(x) and £ (x) must have asymptotic behavior A 2 ~ x ; this is required to obtain a logarithmic dis- 
cretization around the point x — 0. For numerical solution of the equations, it is therefore convenient to use 
the following Ansatz: 

e(x) = g{x)K 2 -\ £(x) = f(x)A 2 ~ x . (34) 
The unknown functions g(x) and f(x) are then 0(1) for all x. The differential equations to be solved are 

^=lnA{g(x)-C[x,g(x)A 2 -*]}, 

fff (a;)A 2 - x , s , (35) 

dx A 2 -*p[f(x)A 2 -*] ' 
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with initial conditions g(2) = 1 and /(l) = 1/A. One first solves for g(x), which is then used in the equation 
for f(x). Both equations are stiff, thus some care is needed in their numeric solution to prevent that g{x) 
or f(x) diverge. 

When NRG is used as an impurity solver in the dynamical mean-field theory, the input to the calculation 
is the effective impurity hybridisation function which contains information on the density of states of the 
sclf-consistently defined medium. The hybridisation function usually takes the form of a tabulated function. 
A differential equations solver was implemented which solves for g(x) and f{x) on a mesh of values of x, 
given an arbitrary input density of states (or hybridisation function). The software package is available from 
the author's home page (http : //nrgljubljana. ijs . si/ adapt). For better portability the solver is written 
in pure ISO CH — h without making use of any external libraries. 

In the solver, p(u>) is calculated for an arbitrary point oj from the tabulated input function by linear 
interpolation. Integrations are performed using trapezoidal method; where the integration boundary falls 
inside a tabulation interval the contribution is, however, calculated by explicitly integrating the linear 
interpolation function. For numerical stability, it is very important that interpolation and integration be 
of consistent order of approximation. The differential equations are solved using fourth-order Runge-Kutta 
solver with adaptive reduction of the integration step size near the boundaries of the tabulation intervals. 
This is especially important when solving for f{x) in situations where p(to) has sharp features (steps, kinks, 
sharp peaks). In the reference implementation of the solver, we use C(x, e) = A/p(e). Constant A is 
first estimated by the numerical integral W = J Q p(u>)dw, but it is then refined by repeatedly solving the 
differential equation for g{x) until g(x) no longer tends to diverge; to determine a suitable value A in this 
shooting-method approach we use the secant method to solve the equation 

g(x max ) = -^r. (36) 

When solving the differential equation for /(&), it is observed that the integration steps need to be made 
small enough in regions of x which correspond to varying p(u>), otherwise f(x) will diverge. For large x in 
the region where p[f(x)A 2 ~ x ] is essentially constant, longer steps may be used. 

The output from the program are tabulated functions f{x) and g(x) which can then be read as input to 
the NRG code. 

6. Conclusion 

A possible improvement would consist of choosing such function C[x, e) to make the discretization grid 
denser in the regions of large variation of p(e), in other words to make it also depend on the derivative 
dp/de. Alternatively, one could systematically study the influence of various environmental modes (i.e. sets 
of nearby conduction-band states) on system dynamics and choose C(x, e) so that sampling is denser in 
energy regions which contribute more [39J. Furthermore, the possibility of making the function C depend 
on the grid parameter x could be used to deform the mesh so that the first weight w\ would not be small 
for small z; w\ tends to with z — > even if the density of states near the band edges is large, which leads 
to similar problems as discussed in Sec. |31 although with lesser severity. 

The author acknowledges discussions with Thomas Pruschke, computer support by GWDG and support 
by the German Science Foundation through SFB 602. 
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